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Abstract 

In dissipative ordinary differential equation systems different time scales cause anisotropic phase volume contraction 
along solution trajectories. Model reduction methods exploit this for simplifying chemical kinetics via a time scale 
separation into fast and slow modes. The aim is to approximate the system dynamics with a dimension-reduced model 
after eliminating the fast modes by enslaving them to the slow ones via computation of a slow attracting manifold. 
We present a novel method for computing approximations of such manifolds using trajectory-based optimization. 
We discuss Riemannian geometry concepts as a basis for suitable optimization criteria characterizing trajectories near 
slow attracting manifolds and thus provide insight into fundamental geometric properties of multiple time scale chem- 
ical kinetics. The optimization criteria correspond to a suitable mathematical formulation of "minimal relaxation" of 
chemical forces along reaction trajectories under given constraints. We present various geometrically motivated cri- 
teria and the results of their application to three test case reaction mechanisms serving as examples. We demonstrate 
that accurate numerical approximations of slow invariant manifolds can be obtained. 

Key words: Model reduction, chemical kinetics, slow invariant manifold, nonlinear optimization, Riemannian 
geometry, curvature 
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1. Introduction 

The need for model reduction in chemical kinetics is mainly motivated by the fact that the computational effort for 
a full simulation of reactive flows, e.g. of fluid transport involving multiple time scale chemical reaction processes, 
is extremely high. For detailed chemical reaction mechanisms involving a large number of chemical species and 
reactions, a simulation in reasonable computing time requires reduced models of chemical kinetics 1 1 1. 

However, model reduction is often also of general interest for theoretical purposes in mathematical modeling. 
Reduced models are intended to describe some essential characteristics of dynamical system behavior while fading 
out other issues. Therefore, they often allow a better insight into complicated reaction pathways, e.g. in biochemical 
systems |2|, and their nonlinear dynamics. 

In dissipative ordinary differential equation systems modeling chemical reaction kinetics different time scales 
cause anisotropic phase volume contraction along solution trajectories. This leads to a bundling of trajectories near 
"manifolds of slow motion" of successively lower dimension as time progresses, illustrated in Figure [T] Many model 
reduction methods exploit this for simplifying chemical kinetics via a time scale separation into fast and slow modes. 
The aim is to approximate the system dynamics with a dimension-reduced model after eliminating the fast modes by 
enslaving them to the slow ones via computation of slow attracting manifolds. 

Very early model reduction approaches like the quasi steady-state (QSSA) and partial equilibrium assumption 
(PEA) 1 1 1 performed "by hand", have set the course for modern numerical model reduction methods that automatically 
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Figure 1: Illustration of trajectories relaxing successively onto a 2-D manifold and a 1-D manifold before converging to equilibrium. Figure 
courtesy of A.N. Al-Khateeb, J.M. Powers, S. Paolucci (private communication). 

compute a reduced model without need for detailed expert knowledge of chemical kinetics by the user. Many of these 
modem techniques are explicitly or implicitly based on a time-scale analysis of the underlying ordinary differential 
equation (ODE) system with the purpose to identify a slow attracting manifold in phase space where - after a short 
initial relaxation period - the system dynamics evolve. For a comprehensive overview see e.g. |3 1 and references 
therein. 

In 1992 Maas and Pope introduced the ILDM-method (H which became very popular and widely used in the 
reactive flow community, in particular in combustion applications. Based on a singular perturbation approach, a local 
time-scale analysis is performed on the Jacobian of the system of ordinary difl'erential equations modeling chemical 
kinetics. Fast time scales are assumed to be fully relaxed and after a suitable coordinate transformation fast variables 
are computed as a function of the slow ones by solving a nonlinear algebraic equation system. For recent developments 
and extensions of the ILDM method see e.g. 1 5 1 and references therein. 

Another popular technique is computational singular perturbation (CSP) method proposed by Lam in 1985 |[6|[3. 
The basic concept of this method is a representation of the dynamical system in a set of "ideal" basis vectors such 
that fast and slow modes are decoupled. In contrast to ILDM and CSP some iterative methods came into application 
that are not based on explicit computation of a time scale separation, but rather an evaluation of functional equations 
suitably describing the central characteristics of a slow attracting manifold, for example invariance and stability. 
Examples are Eraser's algorithm |8, 9, 10 1 and the method of invariant grids EllITKEl. Other widely known and 
successful methods are the constrained runs algorithm |[T3l[T4]|, the rate-controlled constrained equilibrium (RCCE) 
method |15|, the invariant constrained equilibrium edge preimage curve (ICE-PIC) method (I61IT3, and flamelet- 
generated manifolds [IS] [^. In ||2p| finite time Lyapunov exponents and vectors are analyzed for evaluation of 
timescale information. Mitsos et al. formulated an integer linear programming problem explicitely minimizing the 
number of species in the reduced model subject to a given error constraint | 21 1. 

It is obvious that non-local information on phase space dynamics has to be taken into account to get accurate 
approximations of slow attracting manifolds in the general case. Reaction trajectories in phase space that are solutions 
of the ODE system describing chemical kinetics and uniquely determined by their initial values bear such information. 
Based on Lebiedz' idea to search for an extremum principle that distinguishes trajectories on or near slow attracting 
manifolds, we apply an optimization approach for computing such trajectories |[22||23|. Various optimization criteria 
have been suggested |24|, systematically investigated and the trajectory -based approach has been extended to the 
computation of manifolds of arbitrary dimension via parameterized families of trajectories. 

This paper derives and comprehensively discusses various geometrically motivated objective criteria for comput- 
ing trajectories approximating slow attracting manifolds in chemical kinetics as a solution of an optimization problem. 
The corresponding objective functionals are supposed to implicitly incorporate essential characteristics of slow attract- 
ing manifolds related to a minimal remaining relaxation of chemical forces along trajectories on these manifolds. We 
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consider the picture of abstract "dissipative chemical forces" imagined to drive the single elementary reaction steps 
II25II . Due to energy dissipation these forces successively relax while the chemical system is approaching equilibrium. 
The successive relaxation of these forces causes curvature in the reaction trajectories (in the sense of velocity change 
along the trajectory). A slow 1-D manifold in this picture would correspond to a minimally curved reaction trajectory 
along which the remaining relaxation of chemical forces is minimal while approaching chemical equilibrium. 

In particular, we propose and motivate an optimization criterion suitably measuring curvature which is rooted in 
a thermodynamically motivated Riemannian geometry specifically defined for chemical reaction kinetics and based 
on the Second Law of Thermodynamics. This metric provides the phase space of chemical reaction kinetics with a 
geometry specifically capturing the structure of chemical kinetic systems 1 26 J. 

The proposed model reduction method is automatic. The user has to provide only the desired dimension of 
the reduced model and the range of concentrations of the reaction progress variables supposed to parameterize the 
reduced model. For the application examples presented in this work, the numerical optimization algorithm shows fast 
convergence independent of the initial values chosen for numerical initialization. The optimization problems seem to 
be convex for the example systems presented in this paper (see "optimization landscapes" in Section 3.1.2 and 3.3.3 ) 
and would then have a unique solution corresponding to a global minimum of the objective functional. An advantage 
of the presented trajectory optimization approach over local time scale separation methods like QSSA and ILDM is 
the fact that it produces smooth manifolds and whole 1-D manifolds (trajectories) as a solution of a single run of 
the optimization algorithm. Methods based on explicit local time scale separation might yield non-smooth manifolds 
when the fast-slow spectral decomposition changes its structure. Furthermore, the formulation as an optimization 
problem assures results even under conditions where the time scale separation is small and many common model 
reduction methods fail or numerical solutions are difficult to obtain. 



2. Trajectory-based optimization approach 

As described in the introduction, the key of the method presented here is the exploitation of global phase space 
information contained in the behavior of trajectories on their way towards chemical equilibrium. This information 
can be used within an optimization framework for identifying suitable reaction trajectories approximating slow in- 
variant and attracting manifolds (SIM). A suitable formulation as the numerical solution of an optimization problem 
assures the existence of a reduced model irrespective of assumptions on the time scale spectrum, its structure and the 
dimension of the reduced model and sophisticated optimization software can be used for the numerical solution of 
the problem. The central idea behind our approach is that the optimization criterion for the identification of suitable 
trajectories should represent the assumption that chemical forces are - under the given constraints - already maxi- 
mally relaxed along trajectories on the slow attracting manifold. From the opposite point of view this means that the 
remaining relaxation of chemical forces along the trajectories is minimal while approaching chemical equilibrium. 
This means that the velocity change caused by chemical force relaxation is minimal which is intuitively close to the 
notion of a slow manifold. Various ideas for the formulation of suitable optimization criteria are conceivable. 

Mathematically the basic problem can be formulated as 



lin O 

■(0 Jo 

subject to 

dc(0 



min O(c(0) dt (la) 

c(t) 



. =/(c(0) (lb) 
dt 

o = g(cm (ic) 

c,(0) = cl ke /fi,ed. (Id) 

The variables Ck denote the concentrations of chemical species, and /fixed is an index set that contains the indices of 
variables with fixed initial values (the so-called reaction progress variables) chosen to parameterize the reduced model, 
i.e. the slow attracting manifold to be computed. Thus, those variables representing the degrees of freedom within 
the optimization problem are the initial value concentrations of the chemical species q(0), k ^ /fixed- The process of 
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determining c^,k ^ /fixed from c^,k e /fixed is known as "species reconstruction" and represents a function mapping 
the reaction progress variables to the full species composition by determining a point on the slow attracting manifold. 
In our approach, species reconstruction is possible locally, i.e. without being forced to compute the slow attracting 
manifold as a whole. The system dynamics (chemical kinetics determined by the reaction mechanism) are described 
by ([lb]) and enter the optimization problem as equality constraints. Hence an optimal solution of ([T]) always satisfies 
the system dynamics of the full ODE system. Chemical element mass conservation relations that have to be obeyed 
due to the law of mass conservation are collected in the linear function g in ( px| ). The initial concentrations of the 



reaction progress variables are fixed via the equality constraint ([Td]). 

When asymptotically approaching the equilibrium point c^^, which is a stable fixed point attractor in a closed 
chemical system, the system dynamics become infinitely slow and equilibrium will never be reached exactly. By 
approximating the equilibrium point within a surrounding of small radius e > for the concentration of chemical 
species (e.g. the reaction progress variables) by an additional constraint \ck(tf) - c^^\ ^ s (equilibrium composition: 
c^^), the free final time tf can be automatically determined within the optimization problem assuring that this additional 
inequality constraint is fulfilled. However, in practical applications it is usually sufliicient to choose tf large enough 



for the final point of the integration to be close to equilibrium. The objective functional O(c(0) in ( la) characterizes 
the optimization criterion which will be discussed later in detail. 

The key idea of our approach to model reduction in chemical kinetics is found in the fact that suitable trajectories 
can be used to span slow attracting invariant manifolds. The approximated SIM can then be used as a reduced 
model of the underlying ODE model, for example via a look-up table for points on the slow manifold. This reduced 
model is parametrized by the reaction progress variables (coordinate axes) which find a fully natural realization in our 



formulation as trajectory initial concentrations (Id). 



2.1. Numerical Methods: Multiple Shooting in a Parametric Optimization Setting 

The optimization problem ([T]) can be solved as a standard nonlinear optimization problems (NLP), for example 
via the sequential quadratic programming (SQP) method | 27 1. However, one has to decide how to treat the diff'erential 
equation constraint and the objective functional. The easiest way is a decoupled iterative approach, a full numerical 
integration of the ODE model with the current values of the variables subject to optimization. This is called the 
sequential (or single shooting) approach since it fully decouples simulation of the model and optimization. However, 
it is often beneficial to have an "all at once" approach that couples simulation and optimization via discretization of 
the ODE constraint. This simultaneous approach has the advantage of introducing more freedom into the optimization 
problem since the diff'erential equation model does not have to be solved exactly in each iteration of the optimization. 
A beneficial approach to couple the ODE constraint to the optimization is the multiple shooting method. Here, the 
time interval is subdivided into several multiple shooting intervals and additional degrees of freedom are introduced at 
the initial points of each interval. On each multiple shooting interval an independent initial value problem is solved via 
numerical integration. Additional "matching condition" -equality constraints at the level of the optimization problem 
assure continuity of the final solution trajectory between the multiple shooting subintervals. For all results in this 
paper the multiple shooting approach introduced by Bock and Plitt ||28l|29l is used. 

The SQP method basically can be interpreted as Newton's method applied to the Karush-Kuhn-Tucker (KKT) 
necessary optimality conditions of the NLP (see e.g. |30|) and requires the computation of derivatives. For the numer- 
ical approximation of these derivatives by finite diff'erence methods, along with the nominal ODE solution trajectory n 
perturbed trajectories have to be computed, where n is the dimension of the ODE system. To avoid the dependence of 
the resulting derivative on the adaptive discretization schemes of these trajectories provided by an automatic step size 
control in the numerical integration routine, the perturbed trajectories are evaluated on the same grid as the nominal 
trajectory. This approach is called internal numerical diff'erentiation (IND) [31 1. As the systems considered here are 
chemical reaction systems which are usually stiff' systems, the integration itself is performed by DAESOL |32l[^, a 
multistep backward diff'erentiation formula (BDF) diff'erential algebraic equation (DAE) solver. For all computations 
presented in the results section of this paper, the software package MUSCOD-II ||29j [34| |35l has been used. 

For the computation of slow attracting manifolds of dimension larger than one, a sequence of problems of type 
([T]) has to be solved for diff'erent initial values of the reaction progress variables in (Id). For this purpose, we use a 



parametric optimization framework, where neighboring problems are efficiently initialized with the previous optimal 
solution. Through this continuation method embedding the problem into a parametric family of optimization prob- 
lems, the computation of a family of optimal trajectories spanning a higher-dimensional manifold can be signiffcantly 
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accelerated. Such embedding strategy was originally developed and implemented into the package MUSCOD-II by 
Diehl et al. in [ 36 , 31] for fast online optimization, especially real-time optimal control. A variant of this implemen- 
tation has been used for the results presented in this paper. 

2.2. Optimization criteria 

Naturally, the choice of the criterion O(c(0) crucially affects both success and degree of accuracy of the computed 
approximations of the slow attracting manifold. A useful criterion O(c(0) should at least fulfill the following three 
requirements: 

1. O should be physically motivated and describe in a suitable sense the extent of relaxation of "chemical forces" 
or "dynamical modes" in the evolution of reaction trajectories towards equilibrium. 

2. O should be computable from easily accessible data contained in standard models of chemical reaction mecha- 
nisms (e.g. reaction rates, chemical source terms and their derivatives, thermodynamic data). 

3. O should be twice continuously diff'erentiable along reaction trajectories. 

Another desirable but not necessary property, which is related to the invariance of a manifold is the following consis- 
tency property. If a criterion is consistent in the sense of the following definition, the manifold computed as a solution 
of the optimization problem is positively invariant which means that trajectories starting on the manifold at time t^ 
will stay on the manifold for all t to. 

Definition 1 (Consistency property). Suppose an optimal trajectory c(t) has been computed as a solution of ([7]). 
Take the concentrations of the progress variables c(ti) at some time ti > as new initial concentrations for pj| ) and 
solve ([7]) again. Ifc{t) = c(t + t\) holds for the resulting optimal trajectory c(t), we call the optimization criterion O 
consistent. 

The consistency property, illustrated in Fig.|2j can be used as an accuracy test for the computed manifold, because 
the correct attracting SIM should be invariant by definition. However, it poses a strong demand that is not a priori 
incorporated into the problem formulation ([T]) and will not be fulfilled in general for solutions of the optimization 
problem. Nevertheless, an invariant manifold can in principle be constructed in our approach even without a consistent 
criterion by solving M for initial values c^, k e /fixed on the boundary of a desired domain and spanning the low- 
dimensional manifold by the resulting trajectories. 

2.2.7. Curvature -based Relaxation Criteria 

As pointed out before, a suitable optimization criterion O(c(0) should characterize the extent of relaxation of 
"chemical forces". Fundamentally rooted criteria of this type can be derived on the basis of the concept of curvature 
of trajectories in phase space measured in a suitable metric. From a physical point of view curvature is closely related 
to the geometric interpretation of a force and its action on the system dynamics. This picture has a long historical 
tradition. 

One of the most popular examples is Einstein's general theory of relativity |38 | which proposes the idea that 
gravitational force is replaced by a "geometric picture". Einstein's general theory of relativity relates the special 
theory of relativity and Newton's law of universal gravitation with the insight that gravitation can be described by 
curvature of space-time. Space-time is treated as a four-dimensional manifold whose curvature is due to the presence 
of mass, energy, and momentum. 

But even long before Einstein, the concept of curvature has already been related to the concept of force in physics. 
In 1687 Sir Isaac Newton published the laws of motion in his work "Philosophiae Naturalis Principia Mathematica". 
In a diff'erential formulation Newton's second law can be stated sls F = m - a, where m is mass, a is acceleration, and 
F is force. Since the acceleration a is the second derivative of the state variable x(t) with respect to time, a = x, and 
thus contains information about the curvature of x; Newton's law is the first one to directly relate force to curvature. 

In this context it is important to remark that equations of motion in classical mechanics can also be described by a 
variational principle, Hamilton's principle of least action. In Lagrange-Hamilton mechanics 1391 , the trajectory of an 
object is determined in such a way that the action (which is defined as the integral of the Lagrangian over time, where 
the Lagrangian is the diff'erence of kinetic energy and potential energy) is minimal. 
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Figure 2: Illustration of the consistency property: An optimization problem has been solved for a fixed value (=2.0 here) for the progress variable. 
Its solution is the trajectory c(t) starting from the blue circle and converging towards equilibrium (here the coordinate origin, red dot). At a later 
point in time ti > 0, the progress variable is fixed to the value on the trajectory Ck(t) (=1.0 here), k e /fixed, and the optimization problem is solved 
again. If the new solution c(t) coincides with the remaining part of the previous one such as the trajectory starting at the black circle, we call the 
criterion fla) consistent, otherwise (e.g. the trajectory starting at the red circle) it is called inconsistent. 



A central issue in this paper is to transfer the principle of "force = curvature" to the field of chemical systems 
and look for a corresponding variational principle characterizing the kinetics along a slow attracting manifold. In 
chemical systems dissipative forces are active. The diff'erent time scales of dynamic modes result in an anisotropic 
force relaxation in phase space. This force relaxation changes the reaction velocity. Inspired by an analogy observation 
with Newton's geometric interpretation of a force as a second derivative of a trajectory with respect to time, we regard 
the second time-derivative of the chemical composition c(t) characterizing the rate of change of reaction velocity 
through relaxation (dissipation) of chemical forces: 

We consider the tangent (reaction velocity) vectors c(t) = f(c(t)) of reaction trajectories. The relaxation of 
chemical forces results in a change of c(t) along a trajectory on its way towards chemical equilibrium. This change 
along the trajectory may be characterized by taking the directional derivative of the tangent vector of the curve c(t) 
with respect to its own direction v := ^ = ^ . Mathematically this can be formulated as 

D,c(t) := ^mt) + av)\ = Jf(c) • 

da \a=o II/II2 

with Jf(c) being the Jacobian of the right hand side / evaluated at c(t) and || • II2 denoting the Euclidian norm. Hence, 
we may choose the optimization criterion 

\\Jf(c) fh 



in the formulation (la). This criterion bears some resemblance to the recently published method of stretching-based 
diagnostics | 40 | and its application for model reduction (SBR-method). The authors use an expression closely related 
to criterion ^ which measures the stretching of vector fields in the tangent bundle of manifolds. 

The natural way for the evaluation of criterion ^ in the formulation of the objective functional ([la]) would be a 
path integral along the trajectory towards equilibrium 



(D(c(/(0)) d/(0, 

J/(0) 
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where l(t) is the EucHdian length of the curve c(t) at time t given by 
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l(t)= I I|C(T)||2 dr. 

This results in the reparametrization 

d/(0 = ||c(f)ll2 dt. (4) 

The objective used in ( la( would be (using (|2]i): 



rV/(c) /II2 df = f'llclbdr. (5) 

Jo Jo 

However, an alternative norm for the evaluation of \\Jf(c) f\\ might be taken into account, which has already been 
used by Weinhold in [41 1 and is motivated from thermodynamics. This norm is also known as Shashahani norm |42 | 
and is employed for model reduction purposes in 1261 . In this norm the criterion adapted from ^ can be written as 

\\Jf{c)f\\w ^ ifJJ{c)-dmg(l/cd-Jf(c)fy'' 

WfWw (/T.diag(l/c,)-/)i/2 ^' 

with W = diag(l/c/) being the diagonal matrix with diagonal elements 1/c/. This criterion brings thermodynamic 
considerations into play and represents the Riemannian geometry induced by the second differential of Gibbs free 
enthalpy G 

n 

G = Yj ^^[Infe/^') - 1]' ^ = Hess(G). 

It measures the thermodynamic anisotropy of the phase space by weighting the coordinate axis corresponding to 
species / with the gradient ^ = Hess(G)/,/ of the chemical potential jii = of species / in that direction. The 
corresponding metric has been discussed in the context of an entropic scalar product L26J . The corresponding objective 
function in the general optimization problem ([T]) for the W-norm would be 
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\\Jf{c) f\W &t. (7) 



Interestingly, from a different point of view the objective functional ^ and ^ can also be interpreted as mini- 
mizing the length of a trajectory in a suitable Riemannian metric. For any continuously differentiable curve y{t) on a 
Riemannian manifold, the length L of 7 is defined as 



(8) 



with gy(t) being a scalar product defined on the tangent space of the curve in each point. The Riemannian metric gy(t) 
might be chosen as 



8y(t)(f, f) := f Jj(c) ' A . Jf(c) f = \\Jf(c) fWl (9) 



for a positive definite matrix A. The "length-minimizing" objective functional equivalent to criterion ( la) is now 



min L(r) (10) 



subject to constraints ( lb )-( Id). With the solution trajectory of this problem, the "minimum distance from equilibrium 
in a kinetic sense" can be formulated in an explicit mathematical form based on concepts from differential geometry. 
In 1241 . a heuristic choice for a matrix A was made based on the entropy production rate. However, the results 
achieved using the norm proposed in ([6]) yield more accurate results for the computation of slow attracting manifolds 
in chemical kinetics. Another heuristic interpretation of ^ is possible based on the fact that the time-integral over 
a rate of change of velocity is time-averaged velocity whose minimum is intuitively related to the notion of a slow 
manifold. 
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2.2.2. Geometric curvature of curves in space 

In the face of our central aim to relate the principle "force equals curvature" to model reduction for chemical 
kinetics we refer to some mathematical definitions and properties of curvature in this section. Various concepts and 
corresponding definitions of curvature of manifolds and curves in can bee found in general literature on diff'erential 
geometry as e.g. fl3ll44l[45ll . 

Let c : / = (0, L) ^ be a curve defined on an open interval / c R and parameterized by arc length s, meaning 
||£c(5')||2 = 1 (II • lb denoting the Euclidian norm). The value of is a measure for the rate, how rapidly the 

curve pulls away from its tangent in a neighborhood of c{s). That leads directly to the following definition. 

Definition 2 (Curvature). Let c \ J ^ W be a curve parametrized by arc length s e J. The number 



K(s) := 



d^ 



is called curvature of a curve c at s (or at c{s)). 



However, in general trajectories in chemical composition space regarded as curves in vector space are not param- 
eterized in arc-length, but e.g. time t. We want to compute the curvature for the case of an arbitrary parametrization 
t. Let c : / ^ R" be a regular curve, /, / c R open, ip . J ^ I the diff'eomorphism resulting in c := c o ^ being 
parametrized in arc length with w.X.o.g. ip{s) > O'^s e J. Then 



as d(f{s) 



d d 

c((p(s))—(p(s) 
as 



and 



d^ _ _ d^ 
d^^ d(fi(s)^ 



c((f(s))\—(f(s)] -h 



^c(c,(.))-^(.) 



(11) 



(12) 



hold. As c is parametrized in arc length (11) leads to 

d 



d^ 



(p(s) 



For the second derivative of (f, that appears in ([12]), the application of the chain rule yields (with {-rh being the 
Euclidian scalar product) 



d2 iAfifi")) 

ds^""^'^ = 



I d(fi(s) 



Bringing the last two formulae together with ( 12 ) and t = (f(s) we arrive at the formula for the curvature of c(t): 



K(t) := Oc(c(0) := 



-(c(0,c(0)2 



c(t) 



(13) 



11^(011^ i,^vvi,2 

Recalling the discussions in the last section, an alternative optimization criterion for ([la]) could be the curvature 
( [T3] ). In this context the total (integrated) curvature should be the objective function in ([TJ: 



Ktot •= I Ki 

which can be expressed in time-parameterization as 

^(011^(0112 dt 



(s)ds, 



Ktot 



r K(tmt)\\2 dt = r 
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m\\2 



{C{t\c{t))2 



c(t) 



\m\\l 



dt 



(14) 



with K(t) as in equation ( [T3] ). 

Intuitively, the local curvature of a trajectory on its way to equilibrium in phase space should have a peak each 
time it relaxes onto a lower dimensional manifold. Therefore, also this criterion is related to the relaxation of chemical 
forces in some sense. 



2.2.3. Evaluation of the Objective Functional 

From a practical perspective, the computation of the Jacobian for the expression of the different criteria is not 
necessary, as c(t) = Jf(c(t)) f(c(t)) simply is a directional derivative of the ODE vector field with respect to its own 
direction. This directional derivative could be evaluated using classical diff'erence quotients B6ll . but a more appealing 
alternative is found in BTl . Instead of using the central diff'erence formula 

for the approximation of the derivative of the real valued function F(x), Squire and Trapp BTl suggest replacing S 



with i^ (i = V-1). If F is an analytic function, ( 15 ) then reads 



with being the imaginary part of z. This is called complex-step derivative approximation. This result is especially 



appealing, as ( 16 ) does not contain a subtraction and hence eliminates cancellation errors. Therefore d can be chosen 



very small, hence making higher-order terms in the Taylor expansion negligible. For the directional derivative c at a 



point Co with c from c = f(c), reads 

c\co ~ 7 • (17) 



Compared to the use of the full Jacobian, the complexity for the evaluation of c can be reduced from 0(n^) to 0(n) 
using this complex variable approach. At the same time a high accuracy is guaranteed by the possibility of using an 
extremely small 5. 

However, a numerical difficulty occurs for the evaluation of the objectives ^ and ([14]). In case of (|7| the weights 
for the W-norm are obtained as inverted species concentrations. Especially for radical species the denominator be- 
comes generally very small near chemical equilibrium resulting in numerical instabilities. The case of ([14]) is even 
more difficult, as negative exponents > 1 occur for the norm of the reaction rates. Near the equilibrium point the reac- 
tion rates become infinitesimally small and this results in severe numerical problems. A remedy for this problem is an 
additional equality constraint. Instead of fixing the final time tf at a large value, it can be left free in the optimization 
determined by an end point constraint 

ll/(cfe))||2 = € (18) 

with a sufficiently large e keeping the end point of the trajectory away from equilibrium. 
3. Results 

In this section, results for our model reduction method based on trajectory optimization are presented. We choose 
three different chemical reaction mechanisms to demonstrate its application: the Davis-Skodje model system, a sim- 
plified reaction mechanism for the combustion of H2, and a realistic temperature dependent mechanism for ozone 
decomposition. For these three mechanisms all previously discussed objective functionals ([5]), ([7|, and ( [T4] ) in the 
general problem ([TJ are tested for the purpose of numerically computing approximations of slow attracting manifolds. 
The results are compared and discussed. 

3.1. The Davis-Skodje Test Problem 

The well-known Davis-Skodje mechanism is our first test case 1811481. 



dyi dy2 (7 - ^)yi + yyj 

-df = -''^ "dT^-^^^^ (i+,0^ ' 

where 7 > 1 is a measure for the spectral gap or stiffness of the system. Typically model reduction algorithms show a 
good performance for large values of 7, which represent a large gap between the time scales of fast and slow modes. 
Small values of y impose a significantly harder challenge on the computation of slow attracting manifolds. For reasons 
of adjustable time scale separation and analytically computable slow invariant manifold and ILDM, the Davis-Skodje 
model is widely used for testing numerical model reduction approaches. 
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Figure 3: Results for the Davis-Skodje problem with |5j as objective functional. Results for different values of y are shown. The red curve is 
the analytically computed SIM. The black dashed curve represents the analytic Maas-Pope-ILDM. The blue curves are trajectories numerically 
integrated from those initial points that are solutions of our optimization problem. 



We will refer to criterion ^ as criterion A in the following. The dependence of the accuracy of the computed 
slow attracting manifold on the stiffness parameter y becomes obvious. In all cases, the value of yi is fixed as 
reaction progress variable. The optimization problem is solved repeatedly for diff'erent values for yi. For large and 
moderate values of the stiff'ness parameter, good approximations of the SIM (red) are achieved. For y = 1.2 the results 
become more inaccurate. For comparison the Maas-Pope-ILDM is plotted as dashed black line, it can be computed 
analytically for the Davis-Skodje model | 8 |. The c ompu tational eff'ort for one solution of the optimization problem is 
on average about twelve SQP-iterations (cf. Section 2.1 ) which takes about ten seconds in total on a single core Intel® 
Pentium® 4 (3 GHz)-machine with 2 GB memory. Of course, the convergence time (not the result) depends on the 
initial values chosen to start the numerical algorithm. We use a non-equidistant multiple shooting grid with twenty 
intervals. The number of evaluations of the function / is of order 10^, the order of the number of matrix factorizations 
is 10"^. Due to the parametric optimization strategy and application of initial value embedding |36l[37J pointed out in 
Section [2?T| every follow-up solution of a neighboring problem with slightly diff'erent values for the reaction progress 
variables needs only three to five SQP-iterations. 

For the second criterion ^ - denoted B in the following - the results look very similar (see Fig. [4]). As in the 
Davis-Skodje model the values for the variables yi and y2 are of the same order, this is obvious considering the scaling 
in criterion ([6]). Results for criterion C (14), the total curvature, are shown in Figure [5] With this criterion, we could 
not obtain reasonable results for values of y < 6.0 which will be further explained in the next section. 



3.1.2. Optimization Landscapes 

Since the Davis-Skodje test problem consists of only two variables, the structure of optimization landscapes can 
easily be visualized. To compute these landscapes, the initial values of both variables are varied over a fixed range 
and for the trajectories starting in each of these pairs of initial values, the values of (|5]), ([7|, or (14) respectively are 
calculated. These objective functional values are depicted as a function of the initial values of the corresponding 
trajectories. Calculations are performed and plots are generated using MATLAB®. 
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Figure 4: Results for the Davis-Skodje problem with ^ as objective functional. Again, results for different values of y are shown together with 
the SIM (red) and the Maas-Pope-ILDM (black, dashed). 



7 = 6.0 




Figure 5: Results for the Davis-Skodje problem with fT4} as objective function. Only results for y = 6.0 are shown, cf. the discussion in Section 
1 3. 1.2 [ including Figures [6(a)] and |6(b)] 
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For the Davis-Skodje problem, we restrict ourselves to the illustration of criterion C, where we did not achieve 
satisfying results for small values of y (see previous subsection). In Figure 6(a) the results for y = 6.0 are shown. 
Additionally the SIM (red) and the ILDM (black) are projected onto the optimization landscape. Remember, that in 
our approach for a fixed value of yi the minimum of the objective value identifies the corresponding initial value of 

y2- 




Figure 6: Optimization landscape for the Davis-Skodje-model. The integrated (total) curvature fH) is plotted on the z-axis and coded in color for 
illustration reasons. The analytically computed SIM (red) and the analytically given Maas-Pope-ILDM (black) are projected onto the landscape. 



An optimization landscape for y = 2.0 and criterion C is shown in Figure 6(b) In this case, the criterion fails. 



No minimum can be found in the neighborhood of the SIM. The reason for this is that obviously for small time 



scale separation the relation between geometric and kinetic properties of the trajectories pointed out in Section 2.2.2 



becomes weaker and linear segments of trajectories are preferred in the optimization objective functional against a 



slow attracting manifold (see "valley" in Fig. |6(b)| ). As explained in Section [2. 2. 2 [ criterion C measures the curvature 
generated through relaxation of a trajectory onto a slow attracting manifold. Thus, if relaxation is weak, the criterion 
becomes inaccurate. 

3.2. Model Hydrogen Combustion Reaction Mechanism 

In this section we consider a small test mechanism, which has been used for model reduction purposes in ||24l 
[26l 1^ . It consists of six chemical species involved in six (in each case forward and backward) elementary reactions 
involving two element mass conservation relations for hydrogen and oxygen (cf. Table [T]). 

Table 1: Simple hydrogen combustion test mechanism from (26) • Forward and backward rate constants are given temperature-independently. 



Reaction 








k. 


H2 




2H 


2.0 


216.0 


O2 




20 


1.0 


337.5 


H2O 




H + OH 


1.0 


1400.0 


H2 + O 




H + OH 


1000.0 


10800.0 


O2 + H 




+ OH 


1000.0 


33750.0 


H2 + O 




H2O 


100.0 


0.7714 



With the mass conservation relations 



2 + 2 CH20 + c^H + Con = Ci 

2 C02 + CU2O + + <^0H = C2 
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Figure 7: Results for the hydrogen combustion mechanism with criterion A. The free variables and the integrand of the objective functional are 
plotted versus the reaction progress variable ch20- For different values of CH2O the optimization problem is solved and the evolution of the resulting 
trajectories (blue curves) towards equilibrium (red dot) is shown. Especially for the radical species concentrations ch and cq criterion A turns out 
to be inconsistent. 



this mechanism yields a system with four degrees of freedom. For our computations Ci = 2.0 and C2 = 1.0 were 
chosen. 

3.2.1. One -dimensional Manifolds 

We first present results for the computation of one-dimensional manifolds in composition space. The value of 
CH20 serves as reaction progress variable. It is varied between 0.05 and 0.65. We present and compare results for the 



three diff'erent optimization criteria introduced in Section 2.2.1 In general, for testing accuracy of model reduction 



approaches it is very difficult to provide qualitative and quantitative measures of accuracy of numerical results. The 
reason for this is the nonavailability of (analytical) expressions for a SIM. Commonly we consider "eye inspection" 
of trajectory bundling behavior as well as consistency (invariance) as a qualitative measure of accuracy. 

In Figure [7] results for criterion A are shown. The values of the free variables computed in the optimization are 
plotted versus the value of Cu20- Especially for the radical species concentrations ch and cq criterion A turns out to 
be inconsistent to some extent in the sense of Definition [T] Figure [8] shows the results for the weighted criterion B. 
A significant improvement towards better consistency is achieved. The third criterion C which failed in case of the 
Davis-Skodje test problem here performs best, cf. Fi gure [9[ The results are nearly consistent. For criterion C the 
additional equality constraint ([18]), proposed in Section [Z23] has been used to prevent numerical instabilities near the 
equilibrium point. The computational eff'ort for a 1-D manifold is in the order of seconds. 

3.2.2. Two-dimensional Manifolds 

As the hydrogen combustion model has four degrees of freedom, also two-dimensional manifolds can be con- 
structed. In the presented examples, CH20 and ch2 serve as reaction progress variables. We present consistency tests 
plotted in two dimensions and finally show three-dimensional plots of the computed two-dimensional manifold and 
the relaxation of trajectories started from arbitrary initial values onto this 2-D manifold. 



Figure 10 refers to criterion A. The resulting trajectories are plotted. After some time ti the values of the progress 
variables are fixed and the problem is solved again for testing consistency by eye inspection which turns out to be 
quite accurate. For criterion B (see Fig.[TT]), comparably good consistency is observed. 
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Figure 8: Results for the hydrogen combustion mechanism with criterion B. 




0.017 



0.5 1 

CH 2O 



Figure 9: Results for the hydrogen combustion mechanism with criterion C. 
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Figure 10: Results for a two-dimensional manifold for the hydrogen combustion mechanism with criterion A. The values of the progress variables 
are fixed to 0.05 and 0.2 in case of CH2O and to 0.05 in case of • Consistency tests are performed. 




Figure 11: Results for a two-dimensional manifold for the hydrogen combustion mechanism with criterion B. The results are comparable to the 
results for criterion A, cf. Figure [Tol 
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Figure 12: Results of the two-dimensional manifold for the hydrogen combustion mechanism using criterion C. 



Results for the third criterion C are depicted in Figure 12 Here the peaks in the curvature during relaxation onto 
a lower-dimensional manifold mentioned in S ection 17. 2 . 2 1 become obvious. 

For visualization of the two-dimensional manifold, three-dimensional cuts of six-dimensional composition space 



are plotted in Figure [13] The remaining free variables are plotted versus the reaction progress variables. Arbitrary 
trajectories relax on the 2-D manifold spanned by the computed trajectories. The computational effort for a 2-D 
manifold is in the order of some minutes. 



3.3. Ozone Decomposition Reaction Mechanism 

Our last test case is a three component ozone decomposition mechanism (see Table |2]) taken from ISOl . It has 
been chosen to demonstrate the performance of our method taking temperature dependence via Arrhenius kinetics 
into account. Many model reduction approaches explicitly based on time scale separation fail when the spectral gap 
between fast and slow modes becomes too small which is often the case for low temperatures. 



Table 2: Ozone decomposition mechanism from |50|. Rate coefficient k = AT^ Qxp(-Esi/RT). Collision efficiencies in reactions including M: 
fo = 1.14,/o2 = 0.40,/o3 = O^Z 



Reaction 






A 1 cm, mol, s 


b 


E 1^ 


0-hO-hM 




O2 + M 


2.90 X 10^^ 


-1.0 


0.0 


O2 +M 




0-hO-hM 


6.81 X 10^^ 


-1.0 


496.0 


O3 +M 




-h O2 + M 


9.50 X 10^4 


0.0 


95.0 


O2 + M 




O3 -hM 


3.32 X 10^^ 


0.0 


-4.9 






O2 + O2 


5.20 X 10^2 


0.0 


17.4 


O2 + O2 






4.27 X 10^2 


0.0 


413.9 



The ozone decomposition mechanism involves the element mass conservation relation 

Co + 2 + 3 C03 = C 

leaving a system with two degrees of freedom. We choose without loss of generality C = 1. 




Figure 13: Three-dimensional plots of the two-dimensional manifold for the hydrogen combustion mechanism. The free variables are plotted 
versus the progress variables. The manifold is spanned by trajectories (blue) computed for initial fixed values with criterion C as objective 
functional. Arbitrary trajectories (red) fulfilling the element mass conservation are computed to visualize their relaxation onto the manifold. 



3.3.1. Results for Different Optimization Criteria 

The ozone decomposition model has two degrees of freedom and we compute one-dimensional manifolds. The 
results for different temperatures are compared. Consistency tests are performed as in the previous sections. For 
criterion A at T = 1000 K (Fig. 14), a short relaxation phase of the computed trajectories can be observed, indicating 
inconsistency. The other criteria yield more consistent results and seem to be of comparable quality for approximation 
of a slow attracting manifold. 



For a lower temperature ofT = 500 K (Figure 15 ), criterion A obviously fails for the "relatively small" absolute 
values of cq^ , whereas criteria B and C yield good approximations of the SIM. 

For T = 350 K, the effects observed in Figure \T5\ amplify. However, according to the results shown in Figure 16 
criteria B and C still perform reasonably well in this low- temperature region. 



3.3.2. Comparison with ILDM 

For the ozone decomposition mechanism we make a comparison of our results at T 
method |4|. We numerically compute ILDM-points for a range of values. Figure [17] depicts the results. 



1000 K with the ILDM 
A 



comparison of Figure p^b) with the lower row of Figure 14 demonstrates a significantly better performance of our 
trajectory optimization method. The ILDM points do not lie close to the slow attracting manifold. 



3.3.3. Optimization Landscapes 

We compute optimization landscapes for the ozone decomposition model which has two degrees of freedom. 
Initial values for C02 and C03 are varied within the physically allowed range. The value of the objective function is 
computed for trajectories corresponding to tuples of initial values and depicted via color coding in a logarithmic scale. 
We compare these optimization landscapes for T = 1000 K and T = 350 K. 



Figure 18(a) shows the optimization landscape computed for criterion A(T = 1000 K). The other criteria, B and 



C, give rise to a much more distinct minimum of the objective function, cf. Figure 18(b) and 18(c) In the case of 
T = 350 K criterion A fails, cf. Fig. 18(d) no minimum near the SIM is found. The distinct minima for criteria B 
and C become shallow but still allow for an optimal solution close to the SIM. Figures p"8 (e) | and [l 8 (f) | correspond to 



criteria B and C respectively and visualize the results of the optimization problem presented in Section 



3.3.1 
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Co 2 



Co 2 



Co i 



Figure 14: Results for the ozone decomposition mechanism at a temperature ofT= 1000 K with the three different criteria A (left), B (middle), and 
C (right). The free variables are plotted versus cq^ as reaction progress variable. The optimization problem is solved several times with different 
values of co2(0) (x-marks depicting the initial values of the optimal solution trajectories). The solution trajectories starting at the blue x-marks are 
shown on their way to equilibrium (c^^ = 0.5). All criteria work reasonably well, but criterion A is worse concerning consistency. 




Figure 15: Results for the ozone decomposition mechanism at a temperature ofT = 500 K with the three different criteria arranged as in Figure 
[T4] Criteria B and C perform well, whereas criterion A obviously fails. 
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Figure 16: Results for the ozone decomposition mechanism at a temperature ofT = 350 K. 
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(b) 



^0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



Co 2 



Figure 17: Results of the ILDM computation for the ozone decomposition mechanism at a temperature of T = 1000 K. The x-marks depict the 
ILDM-points. They have been computed with the code used in |2 |, accuracy tolerance 10~^ for the solution of the ILDM-equation via Newton's 
method. The computation of ILDM points is initialized with the solution points of the optimization method using criterion B, cf. Fig.[l4] Blue 
lines: trajectories started in the x-marks. 
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(a) Optimization landscape for the ozone decomposition 
mechanism at T = 1000 K. The color represents the value 
of J5j for criterion A. 




(b) Optimization landscape for the ozone mechanism at 
T = 1000 K: The color represents the value of |7j for 
criterion B. 




(c) Optimization landscape for the ozone mechanism at 
T = 1000 K: The color represents the value of fT4} for 
criterion C. 




(d) Optimization landscape for the ozone mechanism at 
T = 350 K: criterion A. 





(e) Optimization landscape for the ozone mechanism at T = 
350 K: criterion B. 



(f) Optimization landscape for the ozone mechanism at T - 
350 K: criterion C. 



Figure 18: Optimization landscape for the ozone decomposition mechanism. The computed value of the different criteria for pairs of trajec- 
tory initial values is depicted in color with the numerical value corresponding to the logarithmically scaled colorbar. The value of cq^ is also 
logarithmically scaled. A trajectory (red) close to the SIM is shown. 
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4. Summary and Discussion 



We present various geometrically motivated criteria for the numerical computation of trajectories approximating 
slow (attracting) invariant manifolds (SIM) in chemical reaction kinetics. The key idea of our approach is to approxi- 
mately span the SIM by trajectories being solutions of an optimization problem for initial values of these trajectories. 
The objective functional is supposed to characterize the extent of relaxation of chemical forces (being minimal in the 
optimal solution) along a reaction trajectory on its way towards equilibrium. Three different criteria are proposed and 
motivated. Whereas the first two criteria use the directional derivative with respect to its own direction of the tangent 
vector field of the kinetic ODE system evaluated in a suitable norm, the third criterion uses a classical diff'erential 
geometric definition of curvature of trajectories regarded as curves in W. 

These criteria are tested with three diff'erent chemical reaction mechanisms: the Davis-Skodje problem |8|, a 
six- species kinetic model for hydrogen combustion, and a realistic ozone decomposition mechanism including tem- 
perature dependence via Arrhenius kinetics. In all cases the quality of the results are evaluated and compared. 

Comparisons with the widely used ILDM-method |4| show that our method bears promise for improvements 
of slow manifold computations in applications. Even though our optimization criteria do not guarantee invariant 
manifolds in general, the solutions in our test examples are close to invariance. It would be possible to compute 
invariant approximations of 1-D manifolds by computing an optimal trajectory for reaction progress variable values 
far from equilibrium und regard the resulting trajectory as a whole as a SIM approximation. Then the manifold would 
be invariant by definition as a trajectory of the ODE system. For the example of a kinetic model for the temperature- 
dependent ozone decomposition it is demonstrated that our approach also works reasonably well in low-temperature 
regions T ^ 1000 K where the ILDM largely fails. 
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